## Require Packages
library(flowCore)
if(!require(flowCore)){
  source("https://bioconductor.org/biocLite.R")
  biocLite("flowCore")
  library(flowCore)
}
library(reshape2)
if(!require(reshape2)){
  install.packages("reshape2")
  library(reshape2)
}
######Ratio CFP and YFP calculation
######Similar to the luciferase assay
ratio_MM <- function(x,c1,c2,log2=TRUE){
  if(log2){
    output <- log2(x[,c1])/log2(x[,c2])
  }else{
    output <- x[,c1]/x[,c2]
  } 
  return(output)
} 
#####Data processing that generates the bin data.
bin.data.FscH<-function(x,bin.size=0.05){
                      file<-read.FCS(x) #Read the fcs into your enviroment
                      file<-exprs(file) #Transform the fcs into a matrix
                      file<-as.data.frame(file) # And into a data.frame
                      colnames(file)<-c("FSC.A", # Giving names to each column 
                                        "FSC.H", # to the data frame from - to dot .
                                        "SSC.A",
                                        "Pacific.Blue", 
                                        "PE.A","Time") 
                      ####Raw data of YFP and CFP these lines select and change the file name
                      Raw.data<-data.frame(file$PE.A,file$Pacific.Blue) 
                      colnames(Raw.data)<-c("YFP","CFP")
                      ####Calculates the ratio and the fold change between CFP and YFP similar
                      #### to luciferase assays
                      Raw.data$ratio<-ratio_MM(Raw.data,2,1,log2=T)
                      Raw.data$fold<-10*(1-Raw.data$ratio)
                      #### Generation of the Analytical function from the raw data using bin size define by bin.:
  
                      ### First the mean of each bin
                      bin<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                   seq(1.5, 10,bin.size)),mean),seq((1.5+bin.size), 10,bin.size))
                      colnames(bin)<-c("CFP","YFP")
                      #### The median of each bin
                      median<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                      seq(1.5, 10,bin.size)),median),seq((1.5+bin.size), 10,bin.size))
                      colnames(median)<-c("median","YFP")
                      #### The max of each bin
                      max<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                   seq(1.5, 10,bin.size)),max),seq((1.5+bin.size), 10,bin.size))
                      colnames(max)<-c("max","YFP")
                      #### The min of each bin
                      min<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                   seq(1.5, 10,bin.size)),min),seq((1.5+bin.size), 10,bin.size))
                      colnames(min)<-c("min","YFP")
                      #### The  Quantile 75
                      qurt75<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                      seq(1.5, 10,bin.size)),quantile, probs=0.75),seq((1.5+bin.size), 10,bin.size))
                      colnames(qurt75)<-c("qurt75","YFP")
                      #### The  Quantile 25
                      qurt25<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                      seq(1.5, 10,bin.size)),quantile, probs=0.25),seq((1.5+bin.size), 10,bin.size))
                      colnames(qurt25)<-c("qurt25","YFP")
                      
                      #### Generating a data frame with all the descriptive statistics 
                      statics<-merge(bin,median,by = "YFP")
                      statics<-merge(statics,max,by="YFP")
                      statics<-merge(statics,min,by="YFP")
                      statics<-merge(statics,qurt75,by="YFP")
                      statics<-merge(statics,qurt25,by="YFP")
                      ###Ratio and fold calculation of binned data 
                      bin$ratio<-(bin$CFP/bin$YFP)
                      bin$fold<-100*(1-bin$ratio)
                      ###All the processed data results together in a list 
                      allinfo<- list(facs=Raw.data, binned= bin,summary=statics)
                      return(allinfo)}
##Same as above but without considering the FSC.A and FSC.H
bin.data<-function(x,bin.size=0.05){file<-read.FCS(x)
                           file<-exprs(file)
                           file<-as.data.frame(file)
                           colnames(file)<-c("FSC.A","SSC.A","Pacific.Blue", "PE.A","Time")
                           Raw.data<-data.frame(file$PE.A,file$Pacific.Blue)
                           colnames(Raw.data)<-c("YFP","CFP")
                           Raw.data$ratio<-ratio_MM(Raw.data,2,1,log2=T)
                           Raw.data$fold<-10*(1-Raw.data$ratio)
                           bin<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                        seq(1.5, 10,bin.size)),mean),seq((1.5+bin.size), 10,bin.size))
                           colnames(bin)<-c("CFP","YFP")
                           median<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                           seq(1.5, 10,bin.size)),median),seq((1.5+bin.size), 10,bin.size))
                           colnames(median)<-c("median","YFP")
                           max<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                        seq(1.5, 10,bin.size)),max),seq((1.5+bin.size), 10,bin.size))
                           colnames(max)<-c("max","YFP")
                           min<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                        seq(1.5, 10,bin.size)),min),seq((1.5+bin.size), 10,bin.size))
                           colnames(min)<-c("min","YFP")
                           qurt75<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                           seq(1.5, 10,bin.size)),quantile, probs=0.75),seq((1.5+bin.size), 10,bin.size))
                           colnames(qurt75)<-c("qurt75","YFP")
                           qurt25<-data.frame(tapply(log2(Raw.data$CFP),cut(log2(Raw.data$YFP),
                                                                           seq(1.5, 10,bin.size)),quantile, probs=0.25),seq((1.5+bin.size), 10,bin.size))
                           colnames(qurt25)<-c("qurt25","YFP")
                           statics<-merge(bin,median,by = "YFP")
                           statics<-merge(statics,max,by="YFP")
                           statics<-merge(statics,min,by="YFP")
                           statics<-merge(statics,qurt75,by="YFP")
                           statics<-merge(statics,qurt25,by="YFP")
                           bin$ratio<-(bin$CFP/bin$YFP)
                           bin$fold<-100*(1-bin$ratio)
                           fit<-lm(formula = bin$CFP ~ bin$YFP+ I(bin$YFP^2))
                           allinfo<- list(facs=Raw.data, binned= bin,fitted=fit,summary=statics)
                           return(allinfo)}
##### Graphic functions: See the main script to understand each part of the function.
plot.miR<-function(data,fil_num =c(1,2,3,4,5),control, 
                   Leg=c("control","1","1","1","1","1"),
                   main= c("titulo"),cex=1,
                   po.pch=19,
                   le.pch=20,
                   mu.cex=0.5,
                   xlim=c(1,10),
                   ylim=c(1,8),
                   l_x=1,
                   l_y=8,
                   col.m=c("black","blue","red","green3","orange","cyan"),
                   xlab="log2 YFP",ylab="log2 CFP"){ 
  data
  plot( data[[control]]$binned$YFP, data[[control]]$binned$CFP, xlim=xlim,ylim=ylim,
        main = main,pch=po.pch,cex=mu.cex,
        xlab=xlab,ylab=ylab,col=col.m[1])
  points( data[[fil_num[1]]]$binned$YFP, data[[fil_num[1]]]$binned$CFP,col=col.m[2],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[2]]]$binned$YFP, data[[fil_num[2]]]$binned$CFP,col=col.m[3],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[3]]]$binned$YFP, data[[fil_num[3]]]$binned$CFP,col=col.m[4],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[4]]]$binned$YFP, data[[fil_num[4]]]$binned$CFP,col=col.m[5],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[5]]]$binned$YFP, data[[fil_num[5]]]$binned$CFP,col=col.m[6],pch=po.pch,cex=mu.cex)
  legend(l_x,l_y, Leg ,col = col.m,cex=cex,pch=le.pch)
  cat("It is so cool to see data plotted, isn't it?")
}
#####The same but plotting the RAtio CFP/YFP
plot.miR.ratio<-function(data,fil_num =c(1,2,3,4,5),control, 
                         Leg=c("control","1","1","1","1","1"),
                         main= c("titulo"),cex=1,
                         po.pch=19,
                         le.pch=20,
                         mu.cex=0.5,
                         xlim=c(1,10),
                         ylim=c(0,1.2),
                         l_x=1,
                         l_y=1.2,
                         col.m=c("black","blue","red","green3","orange","cyan")){ 
  data
  plot( data[[control]]$binned$YFP, data[[control]]$binned$ratio, xlim=xlim,ylim=ylim,
        main = main,pch=po.pch,cex=mu.cex,
        xlab="log2 YFP",ylab="ratio log2(CFP)/log2(YFP)",col=col.m[1])
  points( data[[fil_num[1]]]$binned$YFP, data[[fil_num[1]]]$binned$ratio,col=col.m[2],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[2]]]$binned$YFP, data[[fil_num[2]]]$binned$ratio,col=col.m[3],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[3]]]$binned$YFP, data[[fil_num[3]]]$binned$ratio,col=col.m[4],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[4]]]$binned$YFP, data[[fil_num[4]]]$binned$ratio,col=col.m[5],pch=po.pch,cex=mu.cex)
  points( data[[fil_num[5]]]$binned$YFP, data[[fil_num[5]]]$binned$ratio,col=col.m[6],pch=po.pch,cex=mu.cex)
  legend(l_x,l_y, Leg ,col = col.m,cex=cex,pch=le.pch)
  cat("It is so cool to see data plotted, isn't it?")
}
##Function to export bin data to CSV file to plot in other software
Export<-function(List.FCS,Col.names=fileNames,file="Output.csv"){
  Export<-data.frame(List.FCS[[1]]$binned$YFP) 
  for (i in 1:length(List.FCS)){
    Export[,i+1]<-List.FCS[[i]]$binned$CFP
  }  
  colnames(Export)<-c("YFP",as.character(Col.names))
  write.csv(Export,file)
  return(Export)
}

